/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "dense_perform.hxx"

int main(int argc, char **argv) {

  if (argc != 5) {
    printf("Try :: chol [thread] [blocksize] [ndof] [nitr]\n");
    return 0;
  }

  // ---------------------------------------
  // ** Initialization
  int ndof, blocksize, nthread, nitr;
  nthread   = atoi( (argv[1]) );
  blocksize = atoi( (argv[2]) );
  ndof      = atoi( (argv[3]) );
  nitr      = atoi( (argv[4]) );

  printf("** TEST ENVIRONMENT **\n");
  printf("NDOF      = %d\n", ndof);
  printf("Blocksize = %d\n", blocksize);
  printf("N thread  = %d\n", nthread);
  printf("Iteration = %d\n", nitr);

  int b_mn[2];
  b_mn[0] = b_mn[1] = blocksize;

  FLA_Init();

  // ---------------------------------------
  // ** Environment setting
  double flop = linal::get_flop_chol( 0, ndof );
  FLASH_Queue_set_num_threads(nthread);
  FLASH_Queue_set_sorting(TRUE);
  FLASH_Queue_set_caching(TRUE);

  omp_set_num_threads(nthread);

  // ---------------------------------------
  // ** Matrices
  double 
    t_base, t_flash_repack, t_linal_repack, t_temp, 
    t_flash_decompose, t_linal_decompose;

  linal::Flat_ A;
  A.create(TEST_DATATYPE, ndof, ndof);
  FLA_Random_matrix( ~A );

  FLA_Obj hA_fla;
  t_base = FLA_Clock();
  FLASH_Obj_create_hier_copy_of_flat(~A, 1, (dim_t*)b_mn, &hA_fla);
  t_flash_repack = FLA_Clock()-t_base;

  linal::Hier_ hA_linal;
  t_base = FLA_Clock();
  hA_linal.create(A, blocksize, blocksize);
  t_linal_repack = FLA_Clock()-t_base;

  // ---------------------------------------
  // ** FLASH
  t_flash_decompose = MAX_TIME;
  for (int q=0;q<nitr;++q) {
    printf("*** FLASH::CHOL BEGIN ***\n");
    t_base = FLA_Clock();
    FLASH_Queue_begin();
    FLASH_Chol(FLA_LOWER_TRIANGULAR, hA_fla);
    FLASH_Queue_end();
    t_temp = FLA_Clock()-t_base;
    printf("*** FLASH::CHOL END ***\n");

    t_flash_decompose = min(t_temp, t_flash_decompose);
  }

  // ---------------------------------------
  // ** LINAL
  t_linal_decompose = MAX_TIME;
  for (int q=0;q<nitr;++q) {
    printf("*** LINAL::CHOL BEGIN ***\n");
    t_base = FLA_Clock();
#pragma omp parallel 
    {
#pragma omp single nowait
      linal::dense::chol(FLA_LOWER_TRIANGULAR, hA_linal);
    }
    t_temp = FLA_Clock()-t_base;
    printf("*** LINAL::CHOL END ***\n");

    t_linal_decompose = min(t_temp, t_linal_decompose);
  }

  printf("----------------------------------------------\n");
  printf("*** Report Chol ***\n");
  printf("Ndof       = %d\n", ndof);
  printf("Blocksize  = %d\n", blocksize);
  printf("Nthread    = %d\n", nthread);
  printf("Niteration = %d\n", nitr);
  printf("----------------------------------------------\n");
  printf("Time Supermatrix   = %6.3lf [sec]\n", t_flash_decompose);
  printf("Time LINAL         = %6.3lf [sec]\n", t_linal_decompose);
  printf("----------------------------------------------\n");
  printf("FLOPS Supermatrix  = %6.3lf [Gflops]\n", 
	 flop/t_flash_decompose/1.0e9);
  printf("FLOPS LINAL        = %6.3lf [Gflops]\n", 
	 flop/t_linal_decompose/1.0e9);
  printf("----------------------------------------------\n");
  printf("With repacking cost \n");
  printf("FLOPS Supermatrix  = %6.3lf [Gflops]\n", 
	 flop/(t_flash_decompose+t_flash_repack)/1.0e9);
  printf("FLOPS LINAL        = %6.3lf [Gflops]\n", 
	 flop/(t_linal_decompose+t_linal_repack)/1.0e9);
  printf("----------------------------------------------\n");
  
  // ---------------------------------------
  // ** Matrix
  FLASH_Obj_free(&hA_fla);
  hA_linal.free();
  A.free();

  printf("*** TEST FINISHED ***\n");

  // ---------------------------------------
  // ** Finalization
  FLA_Finalize();
  return 0;
}
